Methods

library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)

df <- read.csv("data/study2.csv") |>
  mutate(
    Date = as.Date(Date),
    Participant = fct_reorder(Participant, Date),
    Screen_Refresh = as.character(Screen_Refresh),
    Illusion_Side = as.factor(Illusion_Side),
    Block = as.factor(Block)
    # Education = fct_relevel(Education, "Master", "Bachelor", "High School", "Other")
  )

Exclusions

# plot(estimate_density(dfraw[dfraw$Participant == "60684f29dbfe1bb2059e5e27_rkqoy", "RT"]))

# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers), which makes it unusable. We understand that you might have been in a hurry or had some other issues; we hope to open-up more slots in the future would you be interested to participate again.

outliers <- c(
  # Error rate of 48.8% Very short RT
  # Prolific Status: NA
  "5e66ed53de7896486f9a71f8_p8ckk"
  )

We removed 1 participants upon inspection of the average error rage (when close to 50%, suggesting random answers) and/or when the reaction time distribution was implausibly fast.

For each block, we computed the error rate and, if more than 50%, we discarded the whole block (as it likely indicates that instructions got mixed up, for instance participants were selecting the smaller instead of the bigger circle).

Error Rate

dfsub <- df |>
  group_by(Participant) |>
  summarize(
    # n = n(),
    Error = sum(Error) / n(),
    RT_Mean = mean(RT),
    RT_SD = sd(RT),
  ) |>
  ungroup() |>
  arrange(desc(Error))

knitr::kable(dfsub) |>
  kableExtra::row_spec(which(dfsub$Participant %in% outliers), background = "#EF9A9A")

Participant	Error	RT_Mean	RT_SD
5e66ed53de7896486f9a71f8_p8ckk	0.488	311	337
613bb8d7741ae6ca04e0431b_gk94a	0.277	921	565
606781fa7584b511b845d52c_cu5nv	0.256	550	152
5d3c6e745602310001bca8aa_6bdtb	0.238	662	495
5c8fcb13dd3b360015b4dff6_dvbq9	0.234	643	327
5e8e06cbdf10df0a811b6f29_g2tm6	0.215	610	293
60127f78468af02bcca0b7a0_lad8d	0.196	576	250
614f8321333460db43b7f229_mpxhl	0.194	1645	983
5ecd58b69c382e0a177e0f61_p3yev	0.162	652	279

Reaction Time Distribution

# RT distribution
p <- estimate_density(df, select = "RT", at = c("Participant", "Block")) |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  mutate(color = ifelse(Participant %in% outliers, "red", "blue")) |>
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = color, group = interaction(Participant, Block), linetype = Block)) +
  geom_vline(xintercept = 150, linetype = "dashed", color = "red") +
  scale_color_manual(values = c("red" = "red", "blue" = "blue"), guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(xlim = c(0, 3000)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  facet_wrap(~Participant) +
  labs(y = "", x = "Reaction Time (ms)")
# p
# ggsave("figures/outliers.png", p, width=20, height=15)

# Filter out
df <- filter(df, !Participant %in% outliers)

Error Rate per Illusion Block

temp <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  summarize(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |>
  arrange(desc(ErrorRate_per_block))

temp2 <- temp |>
  filter(ErrorRate_per_block >= 0.5) |>
  group_by(Illusion_Type, Block) |>
  summarize(n = n()) |>
  arrange(desc(n), Illusion_Type) |>
  ungroup() |>
  mutate(
    n_trials = cumsum(n * 56),
    p_trials = n_trials / nrow(df)
  )

# knitr::kable(temp2)

p1 <- temp |>
  estimate_density(at = c("Illusion_Type", "Block")) |>
  ggplot(aes(x = x, y = y)) +
  geom_line(aes(color = Illusion_Type, linetype = Block)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(y = "Distribution", x = "Error Rate") +
  theme_modern()

p2 <- temp2 |>
  mutate(Block = fct_rev(Block)) |>
  ggplot(aes(x = Illusion_Type, y = p_trials)) +
  geom_bar(stat = "identity", aes(fill = Block)) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
  labs(y = "Percentage of Trials Removed", x = "Illusion Type") +
  theme_modern() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p1 | p2



# Drop
df <- df |>
  group_by(Participant, Illusion_Type, Block) |>
  mutate(ErrorRate_per_block = sum(Error) / n()) |>
  ungroup() |>
  filter(ErrorRate_per_block < 0.5) |>
  select(-ErrorRate_per_block)

rm(temp, temp2)

Reaction Time per Trial

df <- df |> 
  group_by(Participant, Error) |> 
  mutate(Outlier = ifelse(Error == 0 & (RT < 150 | standardize(RT) > 4), TRUE, FALSE)) |> 
  ungroup() 

p1 <- estimate_density(df, select = "RT", at = "Participant") |>
  group_by(Participant) |>
  normalize(select = "y") |>
  ungroup() |>
  merge(df |> 
    filter(Error == 0) |> 
    group_by(Participant) |> 
    summarize(Outlier = mean(RT) + 4 * sd(RT))) |> 
  mutate(Outlier = ifelse(x >= Outlier, TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y)) +
  geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
  geom_line(aes(color = Participant, linetype=Outlier)) +
  geom_vline(xintercept = c(150), linetype = "dashed", color = "red") +
  scale_color_material_d("rainbow", guide = "none") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  guides(linetype = "none") +
  coord_cartesian(xlim = c(0, 4000)) +
  theme_modern() +
  theme(axis.text.y = element_blank()) +
  labs(y = "", x = "Reaction Time (ms)")


p2 <- df |>
  group_by(Participant) |>
  summarize(Outlier = sum(Outlier) / n()) |>
  mutate(Participant = fct_reorder(Participant, Outlier)) |>
  ggplot(aes(x = Participant, y = Outlier)) +
  geom_bar(stat = "identity", aes(fill = Participant)) +
  scale_fill_material_d("rainbow", guide = "none") +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
  see::theme_modern() +
  theme(axis.text.x = element_blank())

p1 | p2

We removed 66 (0.70%) outlier trials (150 ms < RT < 4 SD above mean).

df <- filter(df, Outlier == FALSE)

Participants

dfsub <- df |>
  group_by(Participant) |>
  select(Participant, Age, Sex, Education, Nationality, Ethnicity, Duration, Break_Duration, Screen_Resolution, Screen_Refresh, Device_OS) |>
  slice(1) |>
  ungroup()

The final sample included 8 participants (Mean age = 32.6, SD = 15.2, range: [22, 69]; Sex: 50.0% females, 50.0% males, 0.0% other; Education: Bachelor, 50.00%; High School, 12.50%; Master, 37.50%).

plot_distribution <- function(dfsub, what = "Age", title = what, subtitle = "", fill = "orange") {
  dfsub |>
    ggplot(aes_string(x = what)) +
    geom_density(fill = fill) +
    geom_vline(xintercept = mean(dfsub[[what]]), color = "red", linetype = "dashed") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    ggtitle(title, subtitle = subtitle) +
    theme_modern() +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      plot.subtitle = element_text(face = "italic", hjust = 0.5),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank()
    )
}

plot_waffle <- function(dfsub, what = "Nationality", title = what) {
  ggwaffle::waffle_iron(dfsub, what) |>
    # mutate(label = emojifont::fontawesome('fa-twitter')) |>
    ggplot(aes(x, y, fill = group)) +
    ggwaffle::geom_waffle() +
    # geom_point() +
    # geom_text(aes(label=label), family='fontawesome-webfont', size=4) +
    coord_equal() +
    ggtitle(what) +
    labs(fill = title) +
    theme_void() +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p1 <- plot_distribution(dfsub, "Age", fill = "#FF9800")
p2 <- plot_distribution(dfsub, "Duration", title = "Total Duration", subtitle = "in minutes", fill = "#F44336")
p3 <- plot_distribution(dfsub, "Break_Duration", title = "Break Duration", subtitle = "in minutes", fill = "#3F51B5")

p4 <- plot_waffle(dfsub, "Sex") +
  scale_fill_manual(values = c("Male" = "#2196F3", "Female" = "#E91E63", "Other" = "#FF9800"))

p5 <- plot_waffle(dfsub, "Education") +
  scale_fill_viridis_d()

p6 <- plot_waffle(dfsub, "Nationality") +
  scale_fill_metro_d()

p7 <- plot_waffle(dfsub, "Ethnicity") +
  scale_fill_manual(values = c("Latino" = "#FF5722", "Asian" = "#FF9800", "Caucasian" = "#2196F3", "African" = "#4CAF50", "Jewish" = "#9C27B0"))

p8 <- plot_waffle(dfsub, "Screen_Resolution", title = "Screen Resolution") +
  scale_fill_pizza_d()

p9 <- plot_waffle(dfsub, "Device_OS", title = "Device OS") +
  scale_fill_bluebrown_d()

# p10 <- plot_waffle(dfsub, "Screen_Refresh") +
#   scale_fill_viridis_d()


(p1 / p2 / p3) | (p4 / p5 / p6) | (p7 / p8 / p9)

Models

Delboeuf

Model Selection

best_models <- function(data) {
  models_err <- list()
  models_rt <- list()
  for (k1 in c("", "_log", "_sqrt", "_cbrt")) {
    for (k2 in c("", "_log", "_sqrt", "_cbrt")) {
      for (side in c("", "--SIDE")) {
        name <- paste0("DIFF", k1, "--", "STRENGTH", k2, side)
        message(name)
        f <- paste0(
          "Illusion_Difference",
          k1,
          " * Illusion_Strength",
          k2,
          " + (1|Participant)"
        )

        if (side == "--SIDE") f <- paste0("Illusion_Side * ", f)

        m <- glmmTMB::glmmTMB(as.formula(paste0("Error ~ ", f)),
          data = data, family = "binomial"
        )
        if (performance::check_convergence(m)) {
          models_err[[name]] <- m
        }

        m <- glmmTMB::glmmTMB(as.formula(paste0("RT ~ ", f)),
          data = data
        )
        if (performance::check_convergence(m)) {
          models_rt[[name]] <- m
        }
      }
    }
  }

  # RT
  to_keep <- compare_performance(models_rt, metrics = c("BIC")) |>
    arrange(BIC) |>
    slice(1:5) |>
    pull(Name)

  test <- test_performance(models_rt[to_keep], reference = 1)
  perf <- compare_performance(models_rt[to_keep], metrics = c("BIC", "R2"))
  rt <- merge(perf, test) |>
    arrange(BIC) |>
    select(Name, BIC, R2_marginal, BF) |>
    mutate(
      BF = insight::format_bf(BF, name = ""),
      Model = "RT"
    )

  # Errors
  to_keep <- compare_performance(models_err, metrics = c("BIC")) |>
    arrange(BIC) |>
    slice(1:5) |>
    pull(Name)

  test <- test_performance(models_err[to_keep], reference = 1)
  perf <- compare_performance(models_err[to_keep], metrics = c("BIC", "R2"))

  merge(perf, test) |>
    arrange(BIC) |>
    select(Name, BIC, R2_marginal, BF) |>
    mutate(
      BF = insight::format_bf(BF, name = ""),
      Model = "Errors"
    ) |>
    rbind(rt)
}

data <- filter(df, Illusion_Type == "Delboeuf")

best_models(data)
##                       Name   BIC R2_marginal       BF  Model
## 1           DIFF--STRENGTH   691      0.5218          Errors
## 2       DIFF_log--STRENGTH   691      0.5136  = 0.983 Errors
## 3      DIFF_sqrt--STRENGTH   692      0.5044  = 0.483 Errors
## 4      DIFF_cbrt--STRENGTH   694      0.4986  = 0.204 Errors
## 5     DIFF--STRENGTH--SIDE   695      0.5433  = 0.124 Errors
## 6  DIFF_log--STRENGTH_sqrt 13946      0.0432              RT
## 7   DIFF_log--STRENGTH_log 13946      0.0431  = 0.894     RT
## 8  DIFF_log--STRENGTH_cbrt 13946      0.0431  = 0.886     RT
## 9      DIFF--STRENGTH_sqrt 13947      0.0429  = 0.738     RT
## 10     DIFF--STRENGTH_cbrt 13947      0.0428  = 0.660     RT

Error Rate

Descriptive

plot_descriptive_err <- function(data, side = "leftright") {
  # Sanity checks
  if (length(unique(data$Illusion_Difference)) != 8) stop("Illusion_Difference values != 8")
  if (length(unique(data$Illusion_Strength)) != 15) stop("Illusion_Strength values != 15")

  if (side == "leftright") {
    x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
    x <- tools::toTitleCase(gsub("arrow", "", x))
    if (x == "Left") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Left", "Right")
    } else if (x == "Right") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Right", "Left")
    }
  } else {
    x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
    x <- tools::toTitleCase(gsub("arrow", "", x))
    if (x == "Up") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Up", "Down")
    } else if (x == "Down") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Down", "Up")
    }
    data$Answer <- fct_rev(data$Answer)
  }

  dodge1 <- 0.1 * diff(range(data$Illusion_Difference))
  dodge2 <- -0.1 * diff(range(data$Illusion_Strength))

  p1 <- data |>
    group_by(Illusion_Difference, Illusion_Strength, Answer) |>
    summarize(Error = mean(Error), .groups = "drop") |>
    ungroup() |>
    ggplot(aes(x = Illusion_Difference, y = Error)) +
    geom_bar(aes(fill = Illusion_Strength, group = Illusion_Strength), position = position_dodge(width = dodge1, preserve = "single"), stat = "identity", width=0.03) +
    # geom_line(aes(color = Illusion_Strength, group=Illusion_Strength), position = position_dodge(width=dodge1)) +
    geom_hline(yintercept = 0.5, linetype = "dotted", alpha = 0.3) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_fill_gradientn(colours = c("#4CAF50", "#009688", "#00BCD4", "#2196F3", "#3F51B5", "#673AB7", "#9C27B0")) +
    theme_modern() +
    labs(
      color = "Illusion Strength",
      fill = "Illusion Strength",
      y = "Probability of Error",
      x = "Task Difficulty"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))


  p2 <- data |>
    group_by(Illusion_Difference, Illusion_Strength, Answer) |>
    summarize(Error = mean(Error), .groups = "drop") |>
    ungroup() |>
    # mutate(Illusion_Difference = as.factor(round(Illusion_Difference, 4))) |>
    ggplot(aes(x = Illusion_Strength, y = Error)) +
    geom_hline(yintercept = 0.5, linetype = "dotted", alpha = 0.3) +
    geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.6) +
    geom_bar(aes(fill = Illusion_Difference, group = Illusion_Difference), position = position_dodge(width = dodge2, preserve = "single"), stat = "identity", width=0.1) +
    # geom_line(aes(color = Illusion_Difference), position = position_dodge(width=dodge2)) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_fill_gradientn(colours = c("#F44336", "#FFC107", "#4CAF50")) +
    theme_modern() +
    labs(
      color = "Task Difficulty",
      fill = "Task Difficulty",
      y = "Probability of Error",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))

  if (side == "leftright") {
    p <- ((p1 + facet_wrap(~Answer, ncol = 2, labeller = "label_both")) /
      (p2 + facet_wrap(~Answer, ncol = 2, labeller = "label_both"))) +
      plot_annotation(
        title = paste(data$Illusion_Type[1], "Illusion"),
        theme = theme(plot.title = element_text(face = "bold", hjust = 0.5))
      )
  } else {
    p <- ((p1 + facet_wrap(~Answer, nrow = 2, labeller = "label_both")) |
      (p2 + facet_wrap(~Answer, nrow = 2, labeller = "label_both"))) +
      plot_annotation(
        title = paste(data$Illusion_Type[1], "Illusion"),
        theme = theme(plot.title = element_text(face = "bold", hjust = 0.5))
      )
  }
  p
}


plot_descriptive_err(data)

Model Specification

formula <- brms::bf(
  Error ~ Illusion_Difference * Illusion_Strength + 
    (1 + (Illusion_Difference * Illusion_Strength) | Participant),
  family = "bernoulli"
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_delboeuf_err <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 10.0 seconds.
## Chain 1 finished in 10.5 seconds.
## Chain 2 finished in 10.8 seconds.
## Chain 3 finished in 11.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 10.6 seconds.
## Total execution time: 11.2 seconds.

# parameters::parameters(model)

Model Visualization

plot_model_err <- function(data, model) {
  data <- mutate(data, .dots_side = ifelse(Error == 1, "bottom", "top"))
  
  # Get variables
  vars <- insight::find_predictors(model)$conditional
  vardiff <- vars[1]
  varstrength <- vars[2]
  
  # Get predicted
  pred <- estimate_relation(model,
                            at = vars,
                            length = c(NA, 25))

  # Set colors for lines
  colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data[[vardiff]])))
  diffvals <- as.numeric(as.character(unique(sort(pred[[vardiff]]))))
  names(colors) <- diffvals
  
  # Assign color from the same palette to every observation of data (for geom_dots)
  closest <- diffvals[max.col(-abs(outer(data[[vardiff]], diffvals, "-")))]
  data$color <- colors[as.character(closest)]
  data$color <- fct_reorder(data$color, closest)
  
  # Manual jittering
  xrange <- 0.05*diff(range(data[[varstrength]]))
  data$x <- data[[varstrength]]
  data$x[data$x > 0] <- data$x[data$x > 0] - runif(sum(data$x > 0), 0, xrange)
  data$x[data$x < 0] <- data$x[data$x < 0] + runif(sum(data$x < 0), 0, xrange)
  data$x[round(data$x, 2) == 0] <- data$x[round(data$x, 2) == 0] + runif(sum(round(data$x, 2) == 0), -xrange/2, xrange/2)
  
  
  pred |>
    ggplot(aes_string(x = varstrength, y = "Predicted")) +
    geom_dots(
      data = data,
      aes(x=x, 
          y = Error, 
          group = Error,
          side = .dots_side, 
          # order=as.numeric(color) * ifelse(Error == 1, 1, -1)), 
          order=color), 
      fill = data$color,
      color = NA, 
      alpha=0.5) +
    geom_slab(data=data, aes(y = Error)) +
    geom_ribbon(aes_string(ymin = "CI_low", ymax = "CI_high", fill = vardiff, group = vardiff), alpha = 0.2) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_hline(yintercept = c(0.5), linetype = "dotted", alpha = 0.5) +
    geom_line(aes_string(color = vardiff, group = vardiff)) +
    scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = c("#F44336", "#FFC107", "#4CAF50")) +
    scale_fill_gradientn(colours = c("#F44336", "#FFC107", "#4CAF50")) +
    coord_cartesian(xlim=c(min(data[[varstrength]]), max(data[[varstrength]]))) +
    theme_modern() +
    labs(
      title = paste0(data$Illusion_Type[1], " Illusion"),
      color = "Difficulty", fill = "Difficulty",
      y = "Probability of Error",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p_delboeuf_err <- plot_model_err(data, model_delboeuf_err)
p_delboeuf_err

Reaction Time

Descriptive

plot_descriptive_rt <- function(data, side = "leftright") {
  if (side == "leftright") {
    x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
    x <- tools::toTitleCase(gsub("arrow", "", x))
    if (x == "Left") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Left", "Right")
    } else if (x == "Right") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Right", "Left")
    }
  } else {
    x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
    x <- tools::toTitleCase(gsub("arrow", "", x))
    if (x == "Up") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Up", "Down")
    } else if (x == "Down") {
      data$Answer <- ifelse(data$Illusion_Side == 1, "Down", "Up")
    }
    data$Answer <- fct_rev(data$Answer)
  }

  p1 <- data |>
    group_by(Illusion_Difference, Illusion_Strength, Answer) |>
    summarize(RT = mean(RT), .groups = "drop") |>
    ggplot(aes(x = Illusion_Difference, y = RT)) +
    ggdist::stat_slabinterval(aes(fill = Illusion_Strength, group = Illusion_Strength, color = Illusion_Strength), alpha = 0.5, normalize = "groups") +
    # geom_line(aes(color = Illusion_Strength), position = position_dodge(width=dodge1)) +
    scale_y_continuous(expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = c("#4CAF50", "#009688", "#00BCD4", "#2196F3", "#3F51B5", "#673AB7", "#9C27B0")) +
    scale_fill_gradientn(colours = c("#4CAF50", "#009688", "#00BCD4", "#2196F3", "#3F51B5", "#673AB7", "#9C27B0")) +
    theme_modern() +
    labs(
      color = "Illusion Strength",
      fill = "Illusion Strength",
      y = "Reaction Time (ms)",
      x = "Task Difficulty"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))

  p2 <- data |>
    group_by(Illusion_Difference, Illusion_Strength, Answer) |>
    summarize(RT = mean(RT), .groups = "drop") |>
    ggplot(aes(x = Illusion_Strength, y = RT)) +
    geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.6) +
    ggdist::stat_slabinterval(aes(fill = Illusion_Difference, group = Illusion_Difference, color = Illusion_Difference), alpha = 0.5, normalize = "groups") +
    # geom_line(aes(color = Illusion_Difference), position = position_dodge(width=dodge2)) +
    scale_y_continuous(expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = c("#F44336", "#FFC107", "#4CAF50")) +
    scale_fill_gradientn(colours = c("#F44336", "#FFC107", "#4CAF50")) +
    theme_modern() +
    labs(
      color = "Task Difficulty",
      fill = "Task Difficulty",
      y = "Reaction Time (ms)",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))

  if (side == "leftright") {
    p <- ((p1 + facet_wrap(~Answer, ncol = 2, labeller = "label_both")) /
      (p2 + facet_wrap(~Answer, ncol = 2, labeller = "label_both"))) +
      plot_annotation(
        title = paste(data$Illusion_Type[1], "Illusion"),
        theme = theme(plot.title = element_text(face = "bold", hjust = 0.5))
      )
  } else {
    p <- ((p1 + facet_wrap(~Answer, nrow = 2, labeller = "label_both")) |
      (p2 + facet_wrap(~Answer, nrow = 2, labeller = "label_both"))) +
      plot_annotation(
        title = paste(data$Illusion_Type[1], "Illusion"),
        theme = theme(plot.title = element_text(face = "bold", hjust = 0.5))
      )
  }
  p
}

data <- filter(data, Error == 0)

plot_descriptive_rt(data)

Model Specification

formula <- brms::bf(
  RT ~ Illusion_Difference * Illusion_Strength + 
    (1 + (Illusion_Difference * Illusion_Strength) | Participant)
)

# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)

model_delboeuf_rt <- brms::brm(formula,
  data = data,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 9.7 seconds.
## Chain 2 finished in 10.1 seconds.
## Chain 1 finished in 10.5 seconds.
## Chain 4 finished in 10.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 10.2 seconds.
## Total execution time: 10.5 seconds.

# parameters::parameters(model)

Model Visualization

plot_model_rt <- function(data, model) {
  data <- mutate(data, .dots_side = ifelse(Error == 1, "bottom", "top"))
  
  # Get variables
  vars <- insight::find_predictors(model)$conditional
  vardiff <- vars[1]
  varstrength <- vars[2]
  
  # Get predicted
  pred <- estimate_relation(model,
                            at = vars,
                            length = c(NA, 25))

  # Set colors for lines
  colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data[[vardiff]])))
  diffvals <- as.numeric(as.character(unique(sort(pred[[vardiff]]))))
  names(colors) <- diffvals
  
  # Assign color from the same palette to every observation of data (for geom_dots)
  closest <- diffvals[max.col(-abs(outer(data[[vardiff]], diffvals, "-")))]
  data$color <- colors[as.character(closest)]
  data$color <- fct_reorder(data$color, closest)
  
  p <- pred |>
    ggplot(aes_string(x = varstrength, y = "Predicted")) +
    ggdist::stat_slab(data=data, aes_string(y="RT", group = vardiff, fill=vardiff), alpha=0.5) +
    geom_ribbon(aes_string(ymin = "CI_low", ymax = "CI_high", fill = vardiff, group = vardiff), alpha = 0.2) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_line(aes_string(color = vardiff, group = vardiff)) +
    scale_y_continuous(expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_gradientn(colours = c("#F44336", "#FFC107", "#4CAF50")) +
    scale_fill_gradientn(colours = c("#F44336", "#FFC107", "#4CAF50")) +
    coord_cartesian(xlim=c(min(data[[varstrength]]), max(data[[varstrength]])),
                    ylim=c(150, 2000)) +
    theme_modern() +
    labs(
      title = paste0(data$Illusion_Type[1], " Illusion"),
      color = "Difficulty", fill = "Difficulty",
      y = "Reaction Time (ms)",
      x = "Illusion Strength"
    ) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
  

}

p_delboeuf_rt <- plot_model_rt(data, model_delboeuf_rt)
p_delboeuf_rt

Visualization

plot_all <- function(data, p_err, p_rt) {
  # Get stimuli
  dat <- df |>
    filter(Error == 0) |> 
    filter(Illusion_Type == unique(data$Illusion_Type),
           Answer %in% c("arrowleft", "arrowup")) |>
    select(Stimulus, Illusion_Strength, Illusion_Difference)
  
  dat <- rbind(
    filter(dat, Illusion_Strength == min(Illusion_Strength)) |> 
    filter(Illusion_Difference %in% c(min(Illusion_Difference), max(Illusion_Difference))),
    filter(dat, Illusion_Strength == max(Illusion_Strength)) |> 
    filter(Illusion_Difference %in% c(min(Illusion_Difference), max(Illusion_Difference))),
    filter(dat, Illusion_Difference == min(Illusion_Difference)) |> 
    filter(Illusion_Strength %in% c(min(Illusion_Strength), max(Illusion_Strength))),
    filter(dat, Illusion_Difference == max(Illusion_Difference)) |> 
    filter(Illusion_Strength %in% c(min(Illusion_Strength), max(Illusion_Strength))))
  
  
  img_leftdown <- filter(dat, Illusion_Difference == max(Illusion_Difference)) |> 
    filter(Illusion_Strength == min(Illusion_Strength)) |> 
    pull(Stimulus) |> 
    unique()
  img_rightdown <- filter(dat, Illusion_Difference == max(Illusion_Difference)) |> 
    filter(Illusion_Strength == max(Illusion_Strength)) |> 
    pull(Stimulus) |> 
    unique()
  img_leftup <- filter(dat, Illusion_Strength == min(Illusion_Strength)) |> 
    filter(Illusion_Difference == min(Illusion_Difference)) |> 
    pull(Stimulus) |> 
    unique()
  img_rightup <- filter(dat, Illusion_Strength == max(Illusion_Strength)) |> 
    filter(Illusion_Difference == min(Illusion_Difference)) |> 
    pull(Stimulus) |> 
    unique()
  
  
  img_leftdown <- paste0("study2/stimuli/", img_leftdown, ".png") |> 
    png::readPNG() |> 
    grid::rasterGrob(interpolate=TRUE) |> 
    patchwork::wrap_elements()
  img_rightdown <- paste0("study2/stimuli/", img_rightdown, ".png") |> 
    png::readPNG() |> 
    grid::rasterGrob(interpolate=TRUE) |> 
    patchwork::wrap_elements()
  img_leftup <- paste0("study2/stimuli/", img_leftup, ".png") |> 
    png::readPNG() |>  
    grid::rasterGrob(interpolate=TRUE) |> 
    patchwork::wrap_elements()
  img_rightup <- paste0("study2/stimuli/", img_rightup, ".png") |> 
    png::readPNG() |> 
    grid::rasterGrob(interpolate=TRUE) |> 
    patchwork::wrap_elements()
  
  p <- ((p_err + labs(x = "")) / (p_rt + labs(title = "")) + patchwork::plot_layout(guides="collect"))
  
  # (
  #   (img_leftup / patchwork::plot_spacer() / img_leftdown + patchwork::plot_layout(heights = c(0.5, 2, 0.5))) | 
  #   (patchwork::plot_spacer() / p / patchwork::plot_spacer()  + patchwork::plot_layout(heights = c(0.5, 2, 0.5))) | 
  #   (img_rightup / patchwork::plot_spacer() / img_rightdown  + patchwork::plot_layout(heights = c(0.5, 2, 0.5)))
  #   ) + 
  #   patchwork::plot_layout(widths = c(0.5, 1, 0.5))
  
  (img_leftup | patchwork::plot_spacer() | img_rightup) / p / (img_leftdown | patchwork::plot_spacer() | img_rightdown) + 
    patchwork::plot_layout(heights = c(0.5, 1.5, 0.5))
}

p_delboeuf <- plot_all(data, p_delboeuf_err, p_delboeuf_rt)
p_delboeuf

Ebbinghaus

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Ebbinghaus")

plot_descriptive_err(data)

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Ebbinghaus", Error == 0)

plot_descriptive_rt(data)

Rod and Frame

Model Selection

data <- filter(df, Illusion_Type == "Rod-Frame")

best_models(data)
##                        Name   BIC R2_marginal       BF  Model
## 1  DIFF_sqrt--STRENGTH_sqrt   739      0.4442          Errors
## 2   DIFF_sqrt--STRENGTH_log   740      0.4409  = 0.671 Errors
## 3       DIFF_sqrt--STRENGTH   740      0.4580  = 0.586 Errors
## 4   DIFF_log--STRENGTH_sqrt   741      0.4402  = 0.345 Errors
## 5       DIFF--STRENGTH_sqrt   742      0.4543  = 0.246 Errors
## 6   DIFF_log--STRENGTH_cbrt 14235      0.0299              RT
## 7  DIFF_cbrt--STRENGTH_cbrt 14235      0.0298  = 0.944     RT
## 8    DIFF_log--STRENGTH_log 14235      0.0297  = 0.917     RT
## 9   DIFF_cbrt--STRENGTH_log 14235      0.0296  = 0.853     RT
## 10  DIFF_log--STRENGTH_sqrt 14235      0.0295  = 0.793     RT

Error Rate

Descriptive

plot_descriptive_err(data)

Reaction Time

Descriptive

data <- filter(data, Error == 0)

plot_descriptive_rt(data)

Vertical-Horizontal

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Vertical-Horizontal")

plot_descriptive_err(data)

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Vertical-Horizontal", Error == 0)

plot_descriptive_rt(data)

Zöllner

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Zöllner")

plot_descriptive_err(data)

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Zöllner", Error == 0)

plot_descriptive_rt(data)

White

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "White")

plot_descriptive_err(data)

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "White", Error == 0)

plot_descriptive_rt(data)

Müller-Lyer

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Müller-Lyer")

plot_descriptive_err(data, side = "updown")

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Müller-Lyer", Error == 0)

plot_descriptive_rt(data, side = "updown")

Ponzo

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Ponzo")

plot_descriptive_err(data, side = "updown")

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Ponzo", Error == 0)

plot_descriptive_rt(data, side = "updown")

Poggendorff

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Poggendorff")

plot_descriptive_err(data, side = "updown")

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Poggendorff", Error == 0)

plot_descriptive_rt(data, side = "updown")

Contrast

Error Rate

Descriptive

data <- filter(df, Illusion_Type == "Contrast")

plot_descriptive_err(data, side = "updown")

Reaction Time

Descriptive

data <- filter(df, Illusion_Type == "Contrast", Error == 0)

plot_descriptive_rt(data, side = "updown")

Individual Scores

Delboeuf

Reliability

random <- model_delboeuf_err 
random <- as.data.frame(model_delboeuf_err)
random <- random[str_detect(names(random), "Participant\\[")]
random <- random |> 
  mutate(Draw = 1:nrow(random)) |>
  pivot_longer(-Draw, names_to="Parameter", values_to = "Value") |> 
  separate(Parameter, into = c("Participant", "Parameter"), sep = ",") |> 
  mutate(Participant = str_remove(Participant, "r_Participant\\["),
         Parameter = str_remove(Parameter, "]"))

random |> 
  ggplot(aes(y= Participant, x = Value)) +
  ggdist::stat_slabinterval() +
  facet_wrap(~Parameter)


data <- random |> 
  mutate(Parameter = paste0("Delboeuf_", Parameter)) |> 
  group_by(Participant, Parameter) |> 
  summarize(Score = as.numeric(bayestestR::map_estimate(Value)), .groups="drop") |>
  ungroup() |> 
  pivot_wider(names_from = Parameter, values_from = Score)

Approximation

dat <- data.frame()
for(participant in unique(data$Participant)) {
  freq <- glm(Error ~ Illusion_Difference * Illusion_Strength, 
              family = "binomial",
              data=filter(df, Participant == participant, 
                          Illusion_Type == "Delboeuf"))
  
  parameters::parameters(freq)[1:2] |> 
    as.data.frame() |> 
    mutate(Participant = participant, 
           Parameter = paste0("freq_Delboeuf_", Parameter))
}

Illusion Game Validation (Study 2)

2022-08-06

Introduction

Methods

Exclusions

Error Rate

Reaction Time Distribution

Error Rate per Illusion Block

Reaction Time per Trial

Participants

Models

Delboeuf

Model Selection

Error Rate

Descriptive

Model Specification

Model Visualization

Reaction Time

Descriptive

Model Specification

Model Visualization

Visualization

Ebbinghaus

Error Rate

Descriptive

Reaction Time

Descriptive

Rod and Frame

Model Selection

Error Rate

Descriptive

Reaction Time

Descriptive

Vertical-Horizontal

Error Rate

Descriptive

Reaction Time

Descriptive

Zöllner

Error Rate

Descriptive

Reaction Time

Descriptive

White

Error Rate

Descriptive

Reaction Time

Descriptive

Müller-Lyer

Error Rate

Descriptive

Reaction Time

Descriptive

Ponzo

Error Rate

Descriptive

Reaction Time

Descriptive

Poggendorff

Error Rate

Descriptive

Reaction Time

Descriptive

Contrast

Error Rate

Descriptive

Reaction Time

Descriptive

Individual Scores

Delboeuf

Reliability

Approximation

Structure

Variability